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SUMMARY 
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The multigrid properties of two data reconstruction methods used for achieving second-order 
spatial accuracy when solving the two-dimensional Euler equations are examined. The data recon- 
struction methods are used with an implicit upwind algorithm which uses linearized backward-Euler 
time-differencing. The solution of the resulting linear system is performed by an iterative procedure. 
In the present study only regular quadrilateral grids are considered, so a red-black Gauss-Seidel itera- 
tion is used. Although the Jacobian is approximated by first-order upwind extrapolation, two alterna- 
tive data reconstruction techniques for the flux integral that yield higher-order spatial accuracy at 
steady state are examined. The first method, probably most popular for structured quadrilateral grids, 
is based on estimating the cell gradients using one-dimensional reconstruction along curvilinear coor- 
dinates. The second method is based on Green’s theorem. Analysis and numerical results for the two- 
dimensional Euler equations show that data reconstruction based on Green’s theorem has superior 
multigrid properties as compared to the one-dimensional data reconstruction method. 


INTRODUCTION 


Multigrid methods have become a popular tool for obtaining steady solutions of the Euler or 
Navier-Stokes equations. Although true multigrid performance is difficult to obtain, there is no doubt 
that multigrid methods can significantly decrease the computer time necessary for convergence. 
However, the gain in performance from a single grid algorithm is directly related to the type of 
smoothing operator used on each level. Although explicit methods may be simple to program and 
have a relatively small number of operation counts, the unconditional stability that implicit methods 
offer tends to greatly overcome their disadvantages. In addition, explicit time advancement methods 
generally do not exhibit good smoothing properties when used with higher-order upwind data recon- 
struction techniques for a system of equations. 

In addition to the time advancement technique, the method of flux evaluation plays an important 
role in algorithm efficiency. One commonly used way to achieve higher order accuracy is to recon- 
struct the data on cell faces appropriately using the cell centered data. For grids which consist of log- 
ically rectangular cells, the most popular approach is to use simple one-dimensional curve fitting 
methods such as used by Anderson etal. [1]. The one-dimensional data reconstruction methods have 
been used with great success in two and three-dimensional CFD codes which use grids consisting of 
logically rectangular cells. 


General fluid dynamics problems may require generating grids around complex shapes for which 
it is difficult to generate a single grid consisting of logically rectangular cells. Using multiple-block 
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grids to model complex geometries has been implemented with success using multigrid algorithms 
[2] [3], Another approach for generating grids around complex geometries is to use triangular ele- 
ments. On unstructured triangular grids, however, data reconstruction methods based on Green’s the- 
orem are more prevalent since this does not require interpolation along a coordinate direction. 

In reference [4] the authors presented a single grid stability analysis and numerical experiments 
of several different data reconstruction methods. In this paper, we extend this work to show the effect 
of the data reconstruction on multigrid performance. The Full- Approximation Scheme (FAS) multi- 
grid method has been incorporated into a quadrilateral-based unstructured grid Euler solver using the 
implicit time marching method of reference [5]. 


GOVERNING EQUATIONS 


The governing equations are the time-dependent Euler equations, which express the conservation 
of mass, momentum, and energy for an inviscid gas. The equations are given by 


^ + -<{>FiWQ = 0 

dt AJ t 1 - 

where A is the area of the cell that is bounded by the contour Q with the outward-pointing unit normal 
n. The state vector Q and the flux vectors F are given as 
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where p is the density, u and v are the x and y components of the velocity, e is the energy per unit vol- 
ume, p is the pressure, and U is the velocity in the direction of the outward pointing normal to the cell 


U = h x u + n y v (3) 

The equations are closed with the equation of state for a perfect gas 


p = {y- l )[ e - p( m 2 +v 2 )/2] 
where yis the ratio of specific heats. 


(4) 


TIME ADVANCEMENT ALGORITHM 


The method used for accelerating the solution to steady state is the Full Approximation Scheme 
(FAS) multigrid method. The technique used for smoothing the errors on each grid level is based on 
the scheme described in reference [5] applied to a grid of quadrilateral cells. The method is an 
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implicit upwind algorithm that uses linearized backward-Euler time differencing. The cell-averaged 
solution vector Q is updated at each time level n with the equations 


I/AQ" = -R(Q") 
q" +1 = Q" + AQ" 


(5) 

( 6 ) 


The ooerator R (O") is the discrete approximation to the flux integral in equation (1) at time level In. 

flux^ are evaluated with Van Leer flux-vector splitting [6] and are second-order accurate tf a hn- 
ear data reconstruction method is used. The operator L is written as 


A dR 
V = -£-1 + — - 

At 0Q" 


(7) 


To minimize the bandwidth and maintain block-diagonal dominance of the matrix V , the Jacobian 
3R" /dQ" is approximated by first-order upwind differencing rather than by exactly lmeanz J 

gain the full benefits of an implicit formulation. However, the scope of this w0 ^ is t0 stability 
effects of various data reconstructions to compute the right-hand side of equation (5) J he y 

andsmoothing analysis presented later assumes the linear system is solved exactly at each time step. 

UPWIND STENCILS 


All of the reconstruction stencils used for the right-hand side of equation (5) in this study 
based on MUSCL-type differencing [6]. In this approach, the flux vector F is sp it into two com 

nents 


F n = F(Q\Q-) = F-(Q + ) + F + (<r) 


( 8 ) 


where 


Qface 


= Qaii +e±< Q> 

The values of Q are determined on each side of a cell face by using an 

reconstructing the cell-centered data on each face as shown in figur . P ™«Qidf*red dif 

torn " face values with Van Leer flux-vector splitting [6], The stenc.ls .ha, are constdered drf 

fer in the interpolation operator 0. 

One of the most common methods of data reconstruction for upwind structured flow solvers is to 
interpolate the dan. to the cell face using only the cells along the 

which is perpendicular to the face [1]. Using the cell numbertng shown tn figure 2, a fatntly 
schemes is given by 
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( 10 ) 


where 


Qface — Q 2 + ^ [(1 — K")A_ + (1 + K’)A + ]Q 2 
Qface = Q 3 - ^[(1 + ^)A_ + (1 - KT)A + ]Q 3 


a+Q< = Q,+i - Q, 
A-Q, = Qi - Q,_, 


( 12 ) 

(13) 


These formulas assume the grid has been transformed from physical (x, y ) space to computa- 
tional (£ n) space where the grid spacing (8£ 5 rf) is unity. Using this family of schemes as the inter- 
polation operator results in the flux integration in a cell depending on a total of 9 cells for - 1 < k < 1 
as shown in figure 3. 



Figure 1. Data reconstruction for upwind fluxes 
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Figure 2. Cell Numbering for k Methods 


B - cell being updated 




Figure 3. 9-point stencil 
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We can examine the relation between the discrete equations (10) - (13) and equation (8) by 
expanding the terms of the equations in a Taylor’s series. Examining the interpolation of Q along a £ 


coordinate line, a Taylor’s series expansion about cell 2 is written as 


Qface ^2 ^ ^ | 


(A Q 2 A! 


+ 




2 d? 






(14) 


!f K= 0, a central difference across cell 2 is used to calculate the gradient so that 


©■( Q )=A«||| 


iVCh-Qi 





(15) 


For K= -1, the gradient is approximated using only one-sided information 


0“(Q) = A£ 


dQ 




|(Q2-Ql) 


(16) 


Although not considered in this study, if K= 1/3, the first and second derivatives of equation (14) are 
estimated with central differences which yield a spatially third-order accurate steady-state solutton in 

one dimension. 

The other stencil used in constructing the data on the face is based on Green’s theorem. This was 
used for triangular grids by Barth and Jespersen [7] and Frink [8]. This method of data reconstruction 
was also used by Anderson [5] on triangular grids in conjuction with the implicit scheme shown here. 
The interpolation operator is evaluated in physical (x, y ) space and is written as 

0 ± (Q) = (VQ-r) ± (17) 

where VQ is the average gradient in the cell and is evaluated using Green’s theorem. 


dQ 

dx 

dQ 

dy 


I|(Q)Mn 

ft 

i|(Q)n/Q 

ft 


(18) 


To evaluate this numerically, inverse-distance weighting is used to transfer the cell-averaged data to 
the nodes [8]. 


^ Qcell, 

Qnode — — 4 " 

yi 


(19) 
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where r, is the distance from the i-th cell center to the node. This reduces to simple averaging for uni- 
form grids. Next, the trapezoidal rule is used to integrate around the cell. The x-component is given 


(20 

Here, A is the area of the cell, nodej and node 2 define face,-. As is the length of face,-, and is the x- 
component of the outward pointing unit normal. The data on the cell faces is then determined using 
(17) where the position vector, r, is computed from the cell center to the face center. Using Green’s 
theorem and the trapezoidal rule results in a stencil of 21 cells for the flux integration. The complete 

procedure for determining Q values on the cell faces is shown in figure 4. 
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Figure 4. Data Reconstruction Using Green’s Theorem 


3. Extrapolate to 
Cell Faces 


TRUNCATION ERROR 


A truncation error analysis for the 9-point stencils using k= 0 and k= -1 as well as the 21-point 
stencil has been shown in reference [4] and is summarized here for completeness. The truncation 
error of each of the three stencils is examined by considering the semi-discrete approximation to a 
scalar advection equation with non-negative coefficients a and b. 


du du . du . 

A +a * +4 * =0 


( 21 ) 


This linear equation is a simplified model of the two-dimensional Euler equations. 


Leaving the equation continuous in time, the spatial derivatives are approximated by each stencil 
and expanded in a Taylor series about the point being updated. The 9-point stencil with k = - 1 leads 
to the following equation: 
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(22) 


where bx and by are the grid spacing in the x and y directions, respectively. The difference approxi- 
mation is second order in the grid spacing with a dispersive leading truncation error term. The 
approximation is also dissipative, as can be seen from the fourth-derivative term of the truncation 
error. For an advection velocity that is aligned with the grid (a or b = 0), the dissipative term reduces 
to a fourth derivative in the flow direction. 

For the 9-point stencil with K = 0 we get the following equation: 


du 
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dx dy 12 
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(23) 


This equation differs from equation (22) in the magnitude of the coefficients of the dispersive and dis- 
sipative terms. We expect this difference formula to be less dissipative than the fully-upwind stencil. 


A Taylor series expansion of the 21-point node-averaged stencil for the scalar advection equation 
gives the following: 


du 
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(24) 


This equations looks remarkably similar to equation (23), as the coefficients of the dispersive and dis- 
sipative terms are identical. However, the dissipative term of the 21-point stencil contains cross deriv- 
atives and looks similar to a biharmonic term. Note that even for a grid-aligned advection velocity the 
cross-derivative term does not vanish. We expect that this difference stencil, although of the same for- 
mal accuracy as the 9-point stencil, will be more dissipative. 


STABILITY ANALYSIS 


The basic stability properties of the upwind stencils considered here were examined in reference 
[4]. A Von Neumann analysis is used to examine the stability and convergence properties of the 9- 
point K= 0 and K=- 1 stencils and the 21 -point stencil. For each of the stencils, the equations are dis- 
cretized according to equations (5) to (7). The operator L" is obtained by first-order interpolation in 
all cases, and the right-hand side R (Q”) is obtained with the three second-order stencils. 
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Although the Von Neumann analysis is commonly applied to the scalar advection equation, we 
examine the stability of the system obtained by linearizing the Euler equations about a constant state, 
similar to the work in reference [5]. Applying a Fourier transform in space to the solution vector Q" 
gives the equation 

Q" = z n Q 0 exp(i‘0)exp(i» ^ 

where (j> = kx/8x, \ jf = Ky/Sy are the Fourier modes in the x and y directions, respectively, and z is 
the amplification factor. Substitution of this expression into (5) yields the following equation 

L{(z-1)Q 0 } = -R{Q 0 } (26) 

where L and R are the Fourier symbols of the left- and right-hand-side operators for the constant- 
coefficient problem. Equations (25) and (26) lead to a generalized eigenvalue problem for z. By rear- 
ranging terms, we define the amplification matrix 

6 = I-t ft (27) 

and z is an eigenvalue of G. The amplification matrix is 4 x 4 and complex; a necessary condition for 
stability is that the magnitude of the eigenvalues of 6 are less than one for all <j> and y/. We will refer 
to the amplification factor for a given mode as the magnitude of the largest eigenvalue for that mode. 
The matrix G depends upon four parameters: the Mach number; the flow direction; the CFL number, 
defined here as c8t/8x, where c is the speed of sound; and the cell aspect ratio, 8y/8x. 

The eigenvalue problem was solved numerically for a series of Fourier modes 0 and i/rin the 
range [ -k , k\ . Below we show the amplification factors for a Mach number of 0.8, flow aligned with 
the grid in the x-direction, a CFL number of 100, and a cell aspect-ratio of 1. These results are typical 
of the stability properties of the implicit scheme at other Mach numbers. 

Shown in figure 5 are the amplification factors for the 9-point stencil with k = -1 and k = 0 
for a CFL of 100. This CFL number represents the asymptotic behavior for the three stencils consid- 
ered here as shown in reference [4], Note that the fully-upwind scheme (k = - 1 ) has very poor 
damping of the short- wavelength modes. As CFL -> «> the amplification factor of the 0 = ±K mode 
asymptotically approaches 1. Although unconditionally stable, the scheme is a very poor smoother 
for an FAS multigrid scheme using high CFL numbers. On the other hand, the upwind-biased stencil 
(jc = 0) leads to a scheme with excellent smoothing properties. All the Fourier modes are very well 
damped; in particular, the checkerboard and sawtooth modes have an amplification factor that tends 
to 0 with increasing CFL numbers. This scheme appears to be a very good multigrid smoother. 

By using the 21 -point stencil to discretize the steady-state operator we get even better stability 
properties, as is seen in figure 6. All the high-frequency modes are damped extremely well; the ampli- 
fication factor for <p, y/= ±n has an asymptote of 0, making this operator an excellent choice as a mul- 
tigrid smoother. 

Considering the 9-point, K= 0 stencil and the 21 -point stencil in the case where the flow is skew 
to the grid, we get the results shown in figure 7. In both cases the damping of the short wavelengths is 
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essentially unchanged. The damping of the long wavelengths is worse, however, and the deterioration 
is somewhat more noticeable for the 9-point stencil, particularly for the intermediate wavelengths. 
The 21-point stencil retains its excellent stability properties over a larger range of wavelengths. 



Figure 5. Amplification factors for 9-point stencil, Mach = 0.8, a = 0, CFL - 100: k- -1 (left) 

and k= 0 (right) 



Figure 6. Amplification factor for 21 -point stencil, Mach = 0.8, CL = 0, CFL = 100 

Shown in figure 8 are the smoothing factors, defined as the maximum of the amplification factor 
over the range k/2 < | $ , j i//| < tt, and average amplification factors for the 9-point stencil over a 
range of CFL numbers from 1 to 1024 and fcfrom -1 to l.The Mach number and flow angle are 0.8 
and 45 degrees, respectively. These plots clearly show that the k = 0 stencil has the best smoothing 
properties for the 9-point stencil. 

A comparison of the smoothing and amplification factors for the 21 -point and the 9-point, K = 0 
stencils is shown in figures 9 and 10. Shown in figure 9 are the smoothing and average amplification 
factors for flow aligned with the grid. Note that for CFL numbers up to about 16, the smoothing fac- 
tors are identical. The asymptotic smoothing factors are slightly different: 0.524 and 0.563 for the 21- 
point and 9-point stencils, respectively. In contrast to the smoothing factors, the average amplification 
factor is about 50% lower for the 21 -point stencil compared to the 9-point stencil. In figure 10 plots of 
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Figure 7. Flow at 45 degrees to the grid: 9-point stencil, Mach = 0.8, a = 45 degrees , CFL = 100: 

k= 0 (left) and 21 -point stencil (right) 



Figure 8. Smoothing factors and average amplification factors for K methods 

the smoothing factor and average eigenvalues are shown for both stencils for flow at 45 degrees to the 
grid. The average amplification factors are virtually unchanged, but there is some difference in the 
smoothing factors. The asymptotic values of the smoothing factors have deteriorated, increasing to 
0.554 and 0.628 for the 21-point and 9-point stencils, respectively. The 21-point stencil’s smoothing 
factor is less sensitive to the flow angle than that of the 9-point, k = 0 stencil. 

The effect of grid aspect ratio on the 21 -point and 9-point k= 0 stencil is shown in figure 11. 
Note that there is a large degradation in the smoothing properties for the 9-point K= 0 stencil when 
using high aspect ratio cells such as those in a viscous calculation near a solid wall or wake region. 
The 21-point stencil, however, is generally not affected by the cell aspect ratio. This insensitivity of 
the smoothing factor as the flow angle and grid aspect ratio changes means that we expect that it will 
result in more uniform multigrid performance than the 9-point, k = 0 stencil, over a variety of flow 
conditions and grid topologies. 
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Figure 9. Smoothing factors and average amplification factors for the 21 -point stencil 
and the 9-point stencil, K= 0, for a flow angle of 0 degrees 



Figure 10. Smoothing factors and averag 
and the 9-point stencil, K= 



Log^CFL) 

: amplification factors for the 21 -point stencil 
t, for a flow angle of 45 degrees 


EULER RESULTS 


Results for the two-dimensional Euler equations are now presented. Two test cases are used in 
this study. The first case is the subsonic flow in a channel with a 3% sin 2 * bump. This case was chosen 
because the flow is nearly grid aligned in every cell. The channel length is three times the channel 
height and the length of the bump is equal to the channel height. A freestream Mach number of 0.3 is 
used. The grid used in this study consists of 157 points along the wall and 49 points normal to the 
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L.og, 0 ( dy /dx) 

Figure 1 1 . Effect of Grid Aspect Ratio on Smoothing Factor for Flow Aligned with Grid 

wall and is shown in figure 12. The density contours for the converged solution using the 21 -point 
stencil are also shown in figure 12. All of the cases utilize a 3-level V-cycle using 15 subiterations to 
solve the linear system at each level. One smoothing iteration is performed on each level except the 
coarsest grid where 3 smoothing steps are performed. 

Convergence histories for this case using the 9-point stencil with k= -1 are shown in figure 14a. 
As the CFL increases, the convergence rate improves up to a CFL of about 10 after which the conver- 
gence degrades, eventually becoming unstable. As discussed above, when the CFL is increased, high 
frequency error modes approach neutral stability. The analysis, however, assumes the linear system is 
solved exactly at each time step which is generally not the case with only 15 subiterations. Therefore, 
the scheme may require a prohibitive number of subiterations to remain stable at high CFL numbers. 

The convergence histories for the 9-point stencil with k = 0 are shown in figure 14b. Unlike the 
9-point stencil with k— -1, this stencil produces very good convergence rates as the CFL is increased. 
Note that there is little decrease in the spectral radius after a CFL of 100. This is consistent with the 
analysis shown in figure 10. The convergence histories for the 21 -point stencil are shown in figure 




Figure 12. 3% Sin 2 (x) bump grid and contours 
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14c and are very similar to the 9-point, K= 0 stencil. For this test case, in which the flow is aligned 
with the grid, both stencils have very good convergence properties. 

To examine the behavior of the schemes with higher aspect ratio cells and when the flow is not 
aligned with the grid, a second test case is considered which is a NACA 0012 airfoil in a Mach = 0.8 
freestream at 0 degrees angle-of-attack. The calculations were performed on a 65x25 c-grid which is 
shown in figure 13 along with the converged density contours obtained with the 21 -point stencil. All 
cases were run using a 3-level V-cycle and 20 subiterations to solve the linear system. 



Far Field Grid Near Field Density Contours 


Figure 13. NACA 0012, Mach = 0.8, a = 0° grid and contours 


The convergence histories for both the 21 -point stencil and 9-point stencil with K= 0 are shown 
in figure 14d. Only the k= 0 value is used because of the poor convergence properties of the K= -1 
stencil. As shown, the 21-point stencil converges significantly faster than the 9-point K= 0 stencil. In 
particular, note that the number of multigrid cycles to reach a residual of 10 -16 using the 21-point 
stencil is about the same as for the channel flow. By contrast, the 9-point, k = 0 stencil shows a 
marked deterioration in performance compared to the channel flow case. These results are consistent 
with the analysis for flow angularity and cell aspect ratio effect presented above. 


DISCUSSION 


The analysis and computations presented indicate that the choice of data reconstruction for 
upwind methods can have a substantial effect on the multigrid performance for a given time advance- 
ment scheme. In particular, the popular 9-point, k = -1 stencil exhibits very poor multigrid conver- 
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a) Sin 2 * bump, 9-point, k= -1 stencil 





d) Transonic airfoil, 21 -point and 9-point, 
K= 0 stencil 


Figure 14. Residual Histories 



gence for high CFL numbers. The 9-point, K= 0 stencil has much better smoothing properties but still 
has difficulty damping the high frequency waves if the flow is not aligned with the grid. By using an 
interpolation operator based on Green’s theorem, excellent smoothing properties are obtained for 
high CFL numbers regardless of the flow angularity as shown in figures 9 and 10. This has been 
shown through analysis and confirmed through numerical experiments. 
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